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NUMERICAL APPROXIMATION OF ONE PHASE 
QUADRATURE DOMAINS 

MAHMOUDREZA BAZARGANZADEH AND FARID BOZORGNIA 



jrt Abstract. In this work, we present two numerical schemes for a free boundary 

— i problem called one phase quadrature domain. In the first method by applying 

^ the proprieties of given free boundary problem, we derive a method that leads 

to a fast iterative solver. The iteration procedure is adapted in order to work 
in the case when topology changes. The second method is based on shape 
reconstruction to establish an efficient Shape-Quasi-Newton-Method. Various 
numerical experiments confirm the efficiency of the derived numerical methods. 
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1. Introduction. 

In this paper we shall consider general mathematical approaches to solve 
the free boundary problems of type 



(1.1) 



A(u,n) = o, 



(1.2) 



B(u,n) = 0. 



Here A corresponds to a well posed elliptic boundary value problem in an 
unknown domain 0, = {x : u(x) > 0} = {u > 0} and B operates on 
the functions supported at the free boundary V = dQ. It is supposed that 



function u can be solved from equation (1.1) for any given suitable domain 



VL. More precisely, in this paper we consider the following problem: 



(P) 



u > 0, in 



u 



0, 



m 



>■ > 



\n, 



where ^ is a given measure with compact support. Our aim in this work 
is to study systematic and efficient ways to solve Problem (P) numerically. 
The outline of the paper is as follows. 

In section 2, we present some basic facts and mathematical background of 
quadrature domains. In section 3, we investigate one of the applications of 
quadrature domains, Hele Shaw flow. Section 4 is devoted to derive a numer- 
ical scheme which is based on the properties of the free boundary especially 
blow up techniques. In section 5 we construct another numerical scheme 
for Problem (P) based on shape reconstruction formulation. Finally, in last 
section we investigate some numerical examples which show the efficency of 
the numerical algorithms. 



2. Notations and mathematical background of quadrature 
domains. 

Let us review some notations that we use here. By Q we mean an open subset 
of 1^ and L P (Q) the usual Lebesgue space with respect to the Lebesgue 
measure. HL P {Q) denote the subspace of L p (fi) that consists of harmonic 
functions and SL p (£l) for the subspace of L P {Q) that consists of subharmonic 
functions. We also show the characteristic function of SI by xn and G always 
denotes the usual "fundamental solution" for the Laplace operator in R . 
In other words, for x G K \ {0} 



i 



12-JV 



G(x) 



for N > 3, 

for N = 2, 

where con denotes the volume of the unit ball in R . 



N(N-2)oj n I 
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It is immediately verified that if ft is open and bounded then as function 
of x G ft, 

G(x - y) G IZX 1 (ft), for all y G ft c , 
-G(x-y)eSL 1 (n), for all ye ft. 

Definition 2.1. Let /ibea measure with compact support. By a subhar- 
monic quadrature domain we mean an open connected set ft C M. N such 
that supp(fi) C ft and 



(2.1) / hdx> / hd/i, 

Jn J 

holds for all h G SL l (n). We write ft G Q(/i,S"L 1 ) if (10} holds and 
/x(ft) < oo. 

If one consider J^hdx = J hd/j, for all /i G HL l {Q) then ft is a quadrature 
domain and we write ft G Q^^HL 1 ). 

The simplest quadrature domain is a circular disc. Suppose that fi = a5o 
where 5q is a Dirac mass at origin and a > 0. Then 

Q(/x, .HX 1 ) = QGu, SL 1 ) = {£(0; r)}, 

where r > is determined by |I?(0;r)| = a, (see |6j). Generally if ft is a 
bounded domain in R and 



(2.2) / Mx = |ft| h(x ), 

holds for all h G XL 1 (ft), where xq is an arbitrary point, then ft is a ball 
centered at xq. 

In this article we deal only with subharmonic quadrature domain, and 
from now on by a quadrature domain we mean a subharmonic quadrature 
domain. 

Let U 11 , the Newtonian potential of the measure /i which is defined by 



U»(x) := (G * fi)(x) = G(x- y)dfi(y), x G 



l N , 

and it satisfies the Poisson's equation —AU^ = fi in the distribution sense. 
For the sake of simplicity, we shall use U instead of U xn . 

Gustafsson in [6j has showed that ft G Q(/J,, SL 1 ) if and only if 

(u n <U^, inM^, 
1 ' ) \U cl = Ui i 1 iaR N \Q. 

Also if one considers u = U^ — U > 0, then 
(2.4) Au = X n-V in R N . 



Note that from (2.4) one has Au = xn away from supp(fi) and according 



to the results in local regularity of solutions for elliptic PDEs, we obtain 
rf*(M N ), for every 1< p < oo. Also Vu G W{*\ 



u G W?£(R N ), for every 1< p < oo. Also Vu G W^(R N ). By the Sobolev 
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embedding theorem the first derivatives are therefore Holder continuous with 
Holder exponent a < 1. 

Sakai in [TJ] has proved that the definition of quadrature domain is equiv- 
alent to the well-known one-phase free boundary problem in distribution 
scenes. More precisely, from PDE point of view, Q G Q(u, SL ) is equiva- 
lent to 



(2.5) 



An = xn 
u> 0, 
u = 0, 



M, 



m 
in 
in 



» N 



\n. 



Remark 1. Suppose that m denotes the Lebesgue measure. By (2.3) we 
know that u 







|Vit| = in 
du 



on 



dn 



ds 



I \ Q. Now taking integration of (2.4) gives 

Audx = m(O) — fj,(supp(fj,)) . 




This fact is also a consequence of (2.1 ). In other words, we know the volume 
of the solution priori. 

Example 2.2. As an other example of one phase quadrature domain, sup- 
pose that xq £ ~R N and a > 0, M > 1. Let B\ = B\(xq, a) and u = Mxb%- 
Then (2.5) reads 

/ 26 n |An = xn-Mx Bl , in 0, 

1 ' j |u = |V«| = 0, in O c . 

The spherical symmetry of the problem shows that the we have to find 
a radial solution u = u{\x\) for (2.6). Consequently, we suppose that ft = 
B2 = B2{xq, r) for some r > a. To make more easier we consider that u = u\ 
on B\ and u = u 2 on B 2 \ B\. We desire to patch u\ and u 2 on dB\ without 
loosing regularity. Therefore our problem is 



(2.7) Au = 
with the following conditions 

(2.8) (" 



1-M, inBi, 

1, in B 2 \B U 



U2, Vui = Vu 2 , in _E?i, 
|Vu 2 |=0 in (B 2 



By some calculations and considering the fundamental solution for Laplacian 
operator one has 



(2.9) 



u(x) 



(1-M) 



|z-zo| 2 
2N 



+ A U 



inBi 



|x-xq| 2 
2 AT 



+ A 2 |x-x | +^3, iaB 2 \Bi, 

in (£ 2 ) c , 
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where A\, A 2 , A3 are appropriate constants which are computed with respect 



to (|2.8|). Now we obtain 

;i-M)^ + §^, inB t , 

. 0, in (B 2 ) c , 

where r = M^a. 

For N = 2 we can replace \x — x\ 2 ~ N by log|x — Xq\ in (2.10) and we 
derive that 



(2.11) u(x) = < 



with r = M^a. 



(1 _ M )^Eol_ + ^M+ £(log £ - |), in B ls 

^^ + 4log|x-x |-^(logr + i), mB 2 \B 1 , 
0, in (£ 2 ) c , 



Remark 2. Suppose that N = 2, fj, = 5q. Let -B(0, e) be an approximation 
of supp{p) with M = — 2 f° r e small enough. Then one can obtain r = -j=. 



2.1. An estimate of quadrature domain. 

In Problem (P) the domain £1 is part of the solution and in order to generate 
a mesh, one needs to find a domain which contains Q. To do this we find a 
bigger domain such that £1 is embedded in it as follows. 

By r(/x) we mean a positive number corresponding to the positive measure 
fj, such that 

\ B r(n)\ = m{ B r(^)) = d/j,= /i(l w ), 

J-R N 

where m denotes the Lebesgue measure in M . The following theorem is 
due to Sakai, see [111 . 



Theorem 2.3. p3] Let [i be a finite positive measure with support in the 
closed ball Br, R > 0. Then every quadrature domain fi of /x for subhar- 
monic functions satisfies 

(2.12) n C B r{fl)+R . 

Furthermore, ifr(fi) > 2R then 

B r(fj,)-R C ^- 

For instance let \i = g(x)xB 1 , where g is a positive function with M = 
sup Bl g(x), then C ^ v ^ +r 

For more information about one phase quadrature domain see [H [71 13 [131 
LE]. 
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3. An application (Hele Shaw flow). 

One application of Problem (P) appears in Laplacian growth like Hele-Shaw 
flow which comes up in flow's dynamic. Here we describe the Hele-Shaw 
problem briefly. 

Suppose that some incompressible fluid has been confined between two 
parallel plate and we inject more fluid by moderate velocity to it. Therefore 
the fluid between plates begin to occupy more space. We are interested in 
to study the behavior of the boundary of the fill space. 

To be precise, let i/bea positive, finite and non zero measure with com- 
pact support and supp(v) C D where D is an open subset of M. N by a C 1 
boundary. Let po be a super harmonic function such that 

I —Ann = v in.D, 
I Pd = on oD. 

We are looking for a family of regions Dt for t > 0, such that dDt moves 
with the velocity — Vp.o t . This problem was introduced by S. Richardson 

m. 



Definition 3.1. Suppose that I is an interval in K. Let [i = \D +tu,tE I. 
A map t — Y Dt C M is a weak solution of the free boundary problem if the 
function m € H 1 (R N ) defined by 

(3.2) Aut = XD t ~ V, 

satisfies the following conditions: 

• ut > 0, 

• I M l ~ XD t ) dx = 0. 

Last condition guarantee that ut = in WL N \ Dt, see [5]. 



Remark 3. PDE (3.2) together with the above conditions are a special case 
of Problem Q. 

Next theorem states the corresponding quadrature domain of the solution 
of the Hele-Shaw problem. 

Theorem 3.2. [5j Suppose that \i and Dq are as before and T > 0. Then 
there exists a weak solution 

[0,T] 3t->D t (ZR N , 

for Hele Shaw problem which is unique and if ut be the function appearing 



in (3.2), then ut is also unique and 

'■i 
u t 



/ PD T dr. 
Jo 

Moreover, Dt can be chosen to be 

Dt = D U{z:u t (z)>0}. 
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For more information about Hele shaw see [5], [12], [9] and references 
therein. 



4. First numerical method to approximate the solution of 
Problem (P). 

In this part by applying the properties of free boundary problem, we con- 
struct an algorithm that leads us to a fast iterative solver. The level set 
method is next employed to evolve the interface in the direction of the nor- 
mal velocity field. 

Consider Problem dPl) in dimension one. Our motivation for the first 
method is based on the fact that for any x outside of the supp(^) one has 

u (x) = ±v2u. 

To be more precise, one has 

(4.1) Au = 1, in {x : u(x) > 0} \ supp(fx). 



Let Xf be a free boundary point. Multiply (4.1) by v! and integrate over 
[x , Xf] to find that 

-(u') 2 (x) = u(x). 

Let (c, d) be an initial guess for {x : u(x) > 0} which contains the support 
of measure [i. Next we solve the following boundary value problem 

U" = l-/i in (c,d), 

[ ' ' \u'(c) = ^/Mc), u'{d) = —s/Md). 

Then to get the free boundary points, we move the points c, d in the normal 
direction with speeds y/2u(c) and y/2u(d), i.e, 

df = d + \J 2u(d), Cf = c — \J 2u(c), 

where ct and dt are free boundary points. Note that in this case we need 
only one iteration, see section 4.2. 



Remark 4. To have existence for the boundary value problem (4.2) one has 
to choose (c, d) close enough to supp(^). 

4.1. Blow up techniques and the main idea. 



In higher dimensions we shall prove that when we approach to the free 

boundary still the quotient [ u(x '' goes to one. First, we recall some known 

y/2u(x) 

properties and lemmas that have been proved in [11]. We shall use them in 



the proof of Theorem 4.6 The following lemma shows the growth of u away 



from the free boundary F. 
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Lemma 4.1. [H] Let u G Lf^ c (Vi), 0, = {u > 0} 6e a solution of Problem 
(P). If x G r tfien 

sup it < Cr , 

B r (a:o) 

where C = C(N). 

Corollary 4.2. Lei it 6e as in Lemma \4-1\ Then 

u{x) <Cdist(x,dn) 2 . 

Also we need the following Non degeneracy property of the solutions. 

Lemma 4.3. |11] Let u be a solution of given free boundary problem, then 
we have the inequality 

sup u > — -, for any x G T. 

dB r (x ) »^ 

Definition 4.4. (Local solutions) For given R, M > 0, and xo G I\ let 

Pr(xo,M) be the class of C ' solutions u of Problem (P) in Br(xo) such 
that 

|I>u(a;) — Du(y)\ < M\x — y\, foranyx,yGM . 
In the case x = we also set Pr(M) = P R (0, M). 

In the above definition if R = oo then we get solutions in the entire space 
R and grow quadratically at infinity, which are called global solutions. 

If it G Pr(xq,M) and A > 0, then the proper re-scaling of it at xq is 
defined by 

. . _ u(x + Ax) - u(x ) 

u xo,\[ x ) ~ T2 ' 



Note that by using non degeneracy, Lemma 4.3, and quadratic growth prop- 



erties, Lemma 4.1 it can be shown that when A — > then 

Ux ,x ->■ uo in C lt £(R N ) for any < a < 1, 

where uq G C, ' (M ). This Uq is called a 6/ow; up of it with fixed center xo 
and also uq is a global solution, i.e, uq G Poo(M). For more details see [TT] . 

Theorem 4.5. [11] (Blow up with fixed center). Let it G Pr(xq,M) be a 
solution of Problem (P). Suppose that 

u (x) = lim u XQ> xAx), x G R N , 

j-HX> 

for some sequence Xj — >■ as j — >• oo. T/ien ito ^s homogeneous of degree two 
with respect to the origin, i.e. 

ito(Ax) = A ito(x), for any i£l and A > 0. 

In the proof of next theorem we will use the concept of regular points. A 
xo G T is a regular point if every blow up of it at xo is a half plane solution. 
Precisely, there is two category of blowup for a solution of the Problem (P). 
Let ito be a blowup with a fixed center then it has one of the following forms 
(see HH): 
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Polynomial solution: uo(x) = \{x ■ Ax), x € R . Here A is an 
symmetric matrix with Tr(A) = 1. 



n x n 



• Half plane solutions: uq(x) = \{x ■ e)+, x G 



i N where e is a unit 
vector. 

Theorem 4.6. Let xq be a free boundary point and x £ {u > 0} then 

\Vu(x) 



lim sup = lim inf 

x^x x ^ x y/2u(3 



1. 



Proof. By Theorem |4.5| , blowup solutions at fixed point xq S V is a global 
homogeneous solution of degree two. Let u be a homogeneous global solu- 
tion. Then by above discussion, u has the following form 



u o{x) = -{x- e)\, x £ 



A' 



where e is a unit vector. 



Without loss of generality assume that xq = then we know that 



u(rx) 



(*i)5 



inC 



l.Q 



which means 



and consequently, 



From above one can get 

\2 



u(rx) (xi)'_ 



Vit(rx) 
r 



x\e\ 



0, 



0. 



u(rx) 



(rxij 



+ cr , where c is an arbitrary small constant and 



Vu(rx) = (rxi)ei + 0(r a ), a < 1. 
Using above expression for u{rx) and |Vu(rx)| and taking the quotient im- 
plies the limit. □ 



4.2. A mixed boundary value problem and first algorithm. 

Assume (Q,u) be a smooth solution of Problem (P). Our aim is to build a 
sequence (£l k ,u k ) of solutions of an approximate quadrature domain problem 
which converges towards (D,,u). Assume that p k £ dQ k . Let n k be the 
normal outward vector on d£l k . By Taylor formula, one can write 

d 2 
u(Pk + d k n k ) ~ u(p k ) + d k Vu(p k ) • n fc + y n^ • D 2 u(p k ) ■ n k . 



We wish to have u(p k + n k d k ) = 0, so if we put Vu(p k ) ■ n k = —^2u(p k ) 
and use the approximation D 2 u ~ ^ (Aw)/, then one gets 

(4.3) d k = CVMPk), 

where £ = 2 — \/2- It means that if T^ = d£l k then {F k -\- d k - n k } converges 
to T. We note that in dimension one, d k = yj2u{p k ). 
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To construct an algorithm, let Qq be an initial guess of O which contains 
supp(n). Consider the following boundary value problem which has a vital 
role in the numerical scheme 

\Au = l-fi, infto, 

Remark 5. We note that (4.4) is not stable at the points close to the free 
boundary, therefore alternatively we solve the following problem to have 
more efficient and robust scheme, 

(45) <** »~ onffi, 




We desire that 9u k behaves like \/2uk, therefore one is able to choose 

/ 2 \ 1/2 

(4.6) 0=1 2 - ) . 

V sup Uk-xJ 

dCih-1 

To determine an optimal value for 9 we employe the fixed point process and 
iterate 9 as follows. 



First of all, let 9 be as in (4.6) and solve (4.5) to obtain the value of u on 



dQk- Then compute the corresponding 9 by the formula (4.6) and repeat 



the same process again. This scenario will converge to the best choice of 9. 



The existence of (4.5 ) is based on minimization techniques and is a special 



case of the next lemma. 

Lemma 4.7. [3j Assume (3 : R — > R is smooth, with 

< a < P'{z) <b (ze M), 

for constants a,b. Let f G L 2 (U), U C M^ is a bounded, open set with 
smooth boundary. For 

J — Au = /, in U, 

dn 

there exists a unique weak solution 



^+(3(u) = 0, ondU, 



4.3. Level set formulation. 

The level set method was introduced by Osher and Sethian for implicitly 
tracking dynamic surfaces and curves, see |X0(, I16j . The main idea behind 
this method is to embed an interface T, which lies in R into a surface in 
dimension R . We can do this embedding by defining a proper function </> 
such that r is the zero level set of </>, i.e, 

r = dn = {x e R N ; </>(x) = o}. 
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Suppose that T divides M> N into multiple connected components then one 
can recognize the inside of one component from its exterior when the sign 
of cj> changes. 



Regarding to Theorem 2.3 let T be a given rectangle such that Q C 
B r (u)+R C T for appropriate R > 0. To apply the level set method for 
Problem (P), we need be positive in 7~\ £1 and negative in 0. By this way 
the outward normal vector of fi is given by 

Vcj) 



n 



|V0|" 



We note that Problem (P) is stationary and the level set formulation 
requires a time evolution so we define the parameter t and introduce a 
family of boundaries Q(t) for t > as the level sets by 

dn(t) = {x e R N ; <f)(x,t) = 0}, 
for unknown function (j) : T x M + — > R. By chain rule 

<p t + V^(x(t),t)-x'(t) = 0. 

Let F = x'(t) ■ n which means that F is speed in outward normal direction. 
Then the level set equation will be as follows 

U + F|V0|=O, 
| <f)(x, t = 0) is given. 

In this paper we restrict our attention to the case that eft is considered 
as the sign distance function and therefore \V<f>\ = 1. Hence the level set 
equation turns to 

(4.7) -^ + F = in TxR+. 

Now consider the following boundary value problem 
fA«(t) = l-/i, in n(t), 



According to (4.3), we choose the quantity (^2u(t) as the speed which 



decreases in Q(t) \ supp(n) and goes to zero when Q(t) approaches to the 



free boundary. Regarding to (4.7), the displacement of the boundary fi(£) 



can be obtained by considering the following equation : 

(4.9) ^ + cVMt) = o, on an(t). 

Now let T be the rectangle in section 4.2. The extension of the previous 
equation to whole domain 7", is one of the important issue in the level set 
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approach. To do this we solve the problem: 

'A-y(t) = 1, in U(t) \ supp(n), 

Av(t) = 0, in T\n(t), 

v(t) = 0, on d(supp(p))U dT, 



(4.10) 



We now extend equation (4.9) to T by 



(4.11) -Jl + v (t)=0, mT\supp(n). 

For more information on velocity extension see [4"lll0j. 

4.3.1. First algorithm for Problem (P). 

Choose a tolerance, TOL<< 1. 

(1) Set k = 0, choose an initial domain Qq with To = O^Iq such that 

supp(n) c n C B r{fl)+R . 

(2) Compute u^ on f2/% which is the solution of the following elliptic 
boundary value problem 

( ) i Au k = 1 -^, in ^fc, 



(3) Solve ( 4.10[ ) and obtain v. 



(4) Update the level set function <j) from (4.11) to get r^ + i. 

(5) Solve (*) in &k+i and get Uk+i- 

(6) If sup |u/j +1 | < TOL, then stop else set k = k + 1 and go to (2). 

5. Second numerical method to approach to the solution of 
Problem (P) based on shape optimization. 

The shape sensitivity analysis is used to define a velocity field, which allows 
us to update the surface while decreasing a given cost function. The solution 
of an elliptic boundary value problem usually depends highly nonlinearly on 
the geometry of the given domain. Thus the geometry can not be solved 
straightforward from a linear equation. 

In shape optimization approach, we rewrite the free boundary problem 
such that the minimum of some cost functional is attained at the solution 
of free boundary. The solution of Problem (P) minimizes the functional 

(5.1) E(u,n)= / -\Vu\ 2 dx+ / (l-n)udx, 

Jn 2 Jn 

over u £ H 1 ^) where Q = {u > 0}. Note that we get u = on <9S7. 
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In the following we discuss the shape sensitivity analysis for the above 
shape functional related to Problem (P). At first, we briefly recall some 
basic facts related to shape calculus [20j . 

In shape sensitivity we analyze how the solution of a PDE changes when 
the domain is changing with a velocity field. Let x £ M. N , and V(£,x) be a 
velocity field (vector field) defined inD,VG C k (D; R N ), V\ dD = 0. Let t be 
artificial time. Assume that ECD. It is natural to define transformation 
Tt(V)x = X(t,x) with a velocity field V by differential equations 

dX 

— (t, x) = V(i, x), X(0, x) = x, x e £. 

One can see that this transformation is quite close to a perturbation of the 
identity in [20|, [1] , where the transformation was defined by 

T t (V) = I + tV(x). 

For small perturbations these two transformations are close (see [2T] ). The 
image of S C fi under Tj is £j. 

Let J be a domain functional J : £ i — > R . We say that the functional 
has a directional shape derivative to direction V at £ if the limit 

J(E,)-J(E) 
t-+o t 

exists. If further <iJ(£, V) is linear and continuous with respect to V and 
it exists for all directions V, we say that J is shape differentiable at £. 
By Hadamard's structure theorem, dJ(E,V) depends only on the normal 
component of V on the boundary of E, see [22|. [23] . 

We use the notations uq or u(£l) to show the dependence of solution 
of a given PDE with respect to the domain $7. For a function w(E) and 
E G C , k > 1, we define material derivative as a limit 

v(zj: \){x) := lim . 

This limit may exist either in a weak or a strong sense, and the material 
derivative is called a weak or strong material derivative respectively, see [20j . 
The shape derivative of u(E) in the direction V is the element f'(E; V) 
defined by 

u'(E; V) := t?(E; V) - Vu(S) • V(0), 

whenever it exists either in a weak or a strong sense. For simplicity's sake 
we shall utilize v'-^ instead of f'(E; V). 

Shape derivative represents the change of function v with respect to the 
geometry. Equivalently, shape derivative is the variation of the state variable 
with respect to the shape change. 

The following lemmas represent the basic formulas for shape differentia- 
tion of integrals. In the following we assume that ft is bounded. 
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Lemma 5.1. [20J Let f(Qt) £ L 1 (Qt) be shape differentiable and f'(Clt) 6 
L 1 (Oj), t G [0, T] and T > 0. ijf fit is a C 0,1 -domain, then 

(5.2) (± f f(n t ) dx ) = [ f(n)dx+ [ f(n)<v,n>ds. 



5.1. Shape optimization techniques for Problem (P) and second 
algorithm. 

First ingredient is the shape derivative of the function uq,. 

Lemma 5.2. The shape derivative ofun in the normal direction V, is given 
by the function u'q, satisfies 



(5.3) 



Au' n = 



du 
8 



in £1, 

^<V(0),n>, on dn. 



Proof. The minimizer of the functional in (5.1) satisfies the following equa- 
tion 

Aun t = / = 1 - M in fl*. 

By multiplying a test function, <p £ H^Clt), and taking integral one obtains 



/ Vn^ t -V<p dx = - f ( 



(5.4) / Vn^ f • V99 dx = - / ftpdx. 

Taking the derivative of the above equation respect to t and considering 



Lemma 5.1 one can see that u' n satisfies 



/ Vit'n -Vifdx = - f'dx = 0. 
Jo. " J aa. 



That is 



Au' a = 0. 



The boundary condition in (5.3) is verified by equation (3.6), chapter 3 in 
EDI. □ 



Remark 6. Let T = dQ be the free boundary for the solution of Problem 
(P). Then 

Uq = in $7. 

Let us now to analyze the behavior of the energy near the solution. 



Lemma 5.3. Consider the energy functional (5.1) of Problem (P). Then 
the shape derivative of E with respect to V is 



(5.5) 



dE(£,V) = / div{-\\Vu\ 2 Y)dx. 
Jt, 2 
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IF, 



Proof. By Lemma 5.1 



one can see 



dE(£,V) = f (Vu ■ Vu + (1 - u)u) dx + f (J|Vu| 2 + (1 - u)u)V • nds, 

where v! is the shape derivative of u into direction V. Our assumption on 
Problem (P) states that i*L<9£ = 0. Then the shape derivative of E is 

(5.6) dE(T,, V) = / Vu ■ Vu' dx + / (1 - n)u' + / -\Vu\ 2 V • nds. 

JY, Jt, J BY, 2 



According to Green's theorem, the first term of (5.6) is 
Vu • Vu dx = - 



u'Audx+ I u'——ds, 
on 



ny. 



and we get 



dE(E,V) 



, . , . ,du 

uAudx+ I u 7— as 
an 



c)T 



+ / ' {l-u)u'dx + I \\Vu\ 2 Y-nds 

u'fl — u) dx + / u'^—ds 
Jay. on 



f(l-H)u'dx+ -\Vu\ 2 V-nds 

7s JdT, 2 



,du r f 1,_ ,0,, 
tt~ds + / -VurV-nds. 



^ 



d: 



n 



r^ 



As it is the solution of a Dirichlet problem, Lemma 5.2 gives us u' 
■j^ < V, n > on dTi. Hence we have for dE(Y*, V) the expression 



dE(E,V) = -/ 1 |Vu| 2 V-nds, 
■/as 2 

and by Stock's theorem it turns 

dE(£,V) = f div{- l \Vu\ 2 V)dx. 
Jt, 2 



□ 



Corollary 5.4. T/te solution of Problem (P) is a critical point of the energy 
functional E. 

Proof. We choose V • n = -^ on 9S. If £ C O then ^ < so we 
have dE(E, V) < and it means that E is decreasing respect to V and the 
solution of free boundary where Vu = 0, is a critical point of E. □ 
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5.1.1. Second algorithm for Problem (P). 

(1) Set k = 0, choose an initial domain So such that supp(fi) C So and 
set To = 9So- 

(2) Solve Au k = 1 in Sfe \ supp(fi) with Dirichlet boundary condition 
u k = on r fc , 

(3) Compute a normal velocity from (2), i.e. 

Vn= -Vu k -nr fc . 

(4) Stop if ||VuA;||i,2(r) is sufficiently small. 

(5) Given T k , move the free boundary by Quasi- Newton method, i.e, 

In dimension one 

Xk+i = %k -u'(x k ). 
In dimension two 

T fc+ i = T k - V-u(xfc) • I. 

Obtain the new shape Sfc+i with the free boundary T k+ i. 

(6) Set k = k + 1 and go to (2). 

5.2. Alternative viewpoint. 

One can consider another starting point. We try to determine a shape fi 
such that 

dun n r 

— — = 0, on r. 
an 

In order to derive a suitable weak formulation, we multiply the normal 

derivative by a smooth test function ip and integrate over T, i.e. we have 

dun 



L 



'-tp da = 0. 



/r dn 
By Gauss' Theorem together with the Poisson equation for un we have 

/ (ftp + Vu n ■ V(f) dx = 0, V <p € H%(Sl). 
Jn 

In other words, the first optimality condition for E (with respect to v) reads 

dE(u; if, O) := dE(u + £(p,Q) [ £ =o= / (ff + V-uq • V</j) dx = 0, 

Jn 

for all u S Hq(Q). If one consider 

J (ip, 0.) = / (ftp + Vun ■ V(p) dx, 
Jn 

then J(Q,.) is a continuous linear functional on Hq(Q), i.e, it can be in- 
terpreted as an element of H^ 1 ^) and we can define an operator F(Q) = 
J (CI, .) mapping into H^ 1 ^) such that ( |5.7| ) is equivalent to solving 

(5.7) F(U) = ini? -1 ^). 
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Now we can do all similar calculations for the functional J and deduce same 
results. 

6. Numerical examples. 

Example 6.1. Suppose that /j, = (1 + 2x 2 + y 2 )xB where B is the unit 
ball. We apply the first method and compute the corresponding quadrature 
domain for Problem (P). Let B be the initial guess. Figure nl depicts the 
numerical solution after just three iterations and Figure [2] demonstrates the 
quantity of |Vu| on the boundary of solution. This example confirms that 
the first method is a fast iterative solver. 

Example 6.2. In this example the support of the measure \x is a polygon 
which is shown in Figure [3j We employ the second method and obtain the 
corresponding quadrature domain. Set \i = 1.5%p. In Figure [3j the initial 
guess is the circle and this figure shows the solution after first iteration. 
Figure [4] states the result after four iterations. Figure [5] illustrates the norm 
of the gradient on the boundary of the solution in forth iteration. 

Now let /j, = Hxp- Figure [6] shows the solution after first iteration and 
Figure [7] states the final result which is close to a ball. Figure [9] illustrates 
the quantity of '-j= on a cross section line which has been shown in Figure 
[HI This Figure verifies Theorem |4.6| 



Example 6.3. Suppose that \x = i(xfli +^Xb 2 ) is uniformly distributed on 
two circles Bi(x\,l), B2(x2,l) where x\ = (—2,0), X2 = (\/8, 0). According 
to Example |2.2| or Remark [TJ we find that if t = 4 then B\ and B2 touch 
each other at origin tangentially. Consider the second method and let time 
increase to t = 5 and solve 

-. , j - ~t(XB 1 +2xb 2 ), in tt, 

on <9f2, 




to get the corresponding quadrature domain. Figure 10 shows the solution 



at t = 5 and Figure 11 illustrates |Vu| for t = 5. Figure 12 is the solution 



of similar PDE for t = 6. 
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Surface: u Contour: u 
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Figure 1. The solution of Problem (P) by employing the first 
method and considering /i = (l + 2x 2 +y 2 )xB after three iterations. 
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Figure 2. The quantity |Vu| on the boundary of the solution in 
Figure [TJ 
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Surface: u 
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Figure 3. The colored part shows the solution J7i, after first iter- 
ation, where support of /i is the polygon and the initial guess (S7 ) 
is a ball. 
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Figure 4. Final domain after four iterations when \l = 1.5%p and 
where P is the polygon. 
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Figure 5. The value of |Vu| on the boundary of the solution 
after four iterations. 
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Figure 6. The first iteration for /.i = Hxp, where P is the poly- 
gon. Initial guess is a ball with center at origin. 
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Surface: u 
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Figure 7. Final quadrature domain when /i = ll%p and where 
P is the polygon. 




Figure 8. The surface of the solution u and a cross section line. 
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Figure 9. The amount of '-y= on the cross section line in figure fell. 
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Figure 10. The quadrature domain corresponding to the solution 



of (6.11 for t= 5. 
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Figure 11. The quantity of |Vu| on the boundary of the solution 
of deTTh for t = 5. 



Acknowledgments: This problem was suggested by Professor Henrik Shahgho- 
lian. We wish to thank him for plentiful discussions and profitable sugges- 
tions. 



References 



[1] Delfour M.C, Shape derivatives and Differentiabilty of Min Max, in "shape 
optimization and free boundaries, M.C. Delfour and G. Sabidussi (eds)", Kluwer, 
Dordrecht, 1990, pp 35-339. 

[2] Ebenfelt P., Gustafsson B., Khavinson D., Putinar M., Quadrature Domains and 
Applications, a Harold S. Shapiro Anniversary Volume. , Birkhauser, 2005. 

[3] Evans, Lawrence C. Partial differential equations. Graduate Studies in Mathemat- 
ics 19., American Mathematical Society, Providence, RL, 662 pp, 1998. 

[4] Flusher M., Rumpf M., Bernoulli's free-boundary problem, Qualitative theory and 
numerical approximation, J. Reine Angew. Math., Vol. 486, (1997) 165-204. 

[5] Gustafsson B., Applications of variational inequalities to a moving boundary 
problem for Hele Shaw flows, SIAM J. Math. Anal. 16 (1985), 279-300. 

[6] Gustafsson B., On quadrature domains and an inverse problem in potential theory, 
J.Analyse Math. 55 (1990), 172-216. 



21 



MAHMOUDREZA BAZARGANZADEH AND FARID BOZORGNIA 



Surface: u Contour: u 




-6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 



Figure 12. The quadrature domain corresponding to the solution 



of (6.1) fort = 6. 



[7] Gustafsson B., Putinar M., Selected topics on quadrature domains. Physica D 235, 
90-100 (2007). 

[8] Gustafsson B., Shahgholian H., Existence and geometric properties of solutions of 
a free boundary problem in potential theory, J. Reine Angew. Math. 473 (1996), 
137-179. 

[9] Gustafsson B., Vasilev A., Conformal and potential analysis in Hele-Shaw cells, 
ISBN 3-7643-7703-8, Birkhauser Verlag, 2006. 

[10] Osher J., Fedkiw R., Level set methods and dynamic implicit surfaces, Springer, 
2003. 

[11] Petrosyan A., Shahgholian H., Uraltseva N., Regularity of free boundaries in 
obtacle type problem, book, coming up. 

[12] Richardson S., Helc Shaw flows with a free boundary produced by the injection of 
fluid into a narrow channel, Journal of Fluid Mechanics, 1972 - Cambridge Univ 
Press. 

[13] Sakai M., Quadrature Domains, Lect. Notes Math. Vol. 934, Berlin-Heidelberg: 
Springer- Verlag, 1982. 



[14] Sakai M., Sharp estimates of the distance from a fixed point to the frontier of a 
Hele-Shaw flow. Potential Anal. 8 (1998), no. 3, 277-302. 



NUMERICAL APPROXIMATION OF ONE PHASE QUADRATURE DOMAINS 



25 



[15] Sakai M., Applications of variational inequalities to the existence theorem on 
quadrature domains, Trans. Am. Math. soc. 276(1983), 267-279. 

[16] Sethian J. A, Level set methods and fast marching method, Cambridge University 
Press, 378 pages, 1993. 

[17] Shahgholian H., On quadrature domains and the Schwarz potential, J. Math. Anal, 
and Appl. 171 (1992), 61-78. 

[18] Shahgholian H., Convexity and uniqueness in an inverse problem of potential 
theory, Proc. Amer. Math. Soc. 116 (1992), 1097-1100. 

[19] Shahgholian H., Uraltseva N., Weiss C, The two-phase membrane problem- 
Regularity of the free boundaries in higher dimensions. Int. Math. Res. Not. IMRN 
2007, no. 8. 

[20] Sokolowski J., Zolesio J., Introduction to shape optimization: shape sensitivity 
analysis, Springer, 1992. 

[21] Tiihonen T., Shape optimization and trial methods for free boundary problems, 
Mathematical Modelling and Numerical Analysis, 31:805-825, 1997. 

[22] Tiihonen T., Finite clement approximation of non local heat radiation problems, 
Math. Mod. Meth. Appl. Sci, 1998. 

[23] Zolesio J., Identification dc domaines par deformations, These detat, Univ. Nic, 
1979. 



Department of Mathematics, Stockholm University, S-10691, Stockholm, Sweden 
E-mail address: mahmoudrezaamath.su.se 

Department of Mathematics, Instituto Superior Tecnico, Lisbon. 
E-mail address: bozorgSmath.lst.utl.pt 



